from numpy import matrix
from numpy import linalg
import numpy as np

#------------------------------------------------------------------------------------------------
#Useful procedures:

def inversebigmatrix(B, q):
    """Invert upper triangular matrices with 3x3 matrix-entries and identities on the diagonal, modulo q"""
    n=len(B)
    A=matrix ([[0 for i in xrange(n)] for k in xrange(n)],dtype=object)
    Ainv=matrix ([[0 for i in xrange(n)] for k in xrange(n)],dtype=object)
    for i in xrange(n):
        for j in xrange(n):
            A[i,j]=B[i,j]
            if i == j:
                Ainv[i,j]=matrix([[1,0,0],[0,1,0],[0,0,1]])
            else:
                Ainv[i,j]=matrix([[0,0,0],[0,0,0],[0,0,0]])          
        for k in xrange(n):
            if (i != k):
                factor = A[k, i]
                for l in xrange(n):
                    A[k,l] = (A[k,l] - factor * A[i,l]) % q
                    Ainv[k,l] = (Ainv[k,l] - factor * Ainv[i,l]) % q
    return Ainv

def transfmn(M):
     """Transform a 3x3 matrix into a number"""
     M=M%2
     return 256*M.item(0)+128*M.item(1)+64*M.item(2)+32*M.item(3)+16*M.item(4)+8*M.item(5)+4*M.item(6)+2*M.item(7)+M.item(8)
    
def transfnm(n):
    """Transform a number into a 3x3 matrix"""
    m1 = n
    u11,m2 = m1/256,m1%256
    u12,m1 = m2/128,m2%128
    u13,m2 = m1/64,m1%64
    u21,m1 = m2/32,m2%32
    u22,m2 = m1/16,m1%16
    u23,m1 = m2/8,m2%8
    u31,m2 = m1/4,m1%4
    u32 = m2/2
    u33 = m2%2
    return matrix([[u11,u12,u13], [u21,u22,u23], [u31,u32,u33]])

def extract(k,M):
    """Extracts from a 3x9 matrix the kth 3x3 matrix"""
    if 1<=k<=3:
        uk = matrix([[M.item((0,3*(k-1))), M.item((0,3*(k-1)+1)), M.item((0,3*(k-1)+2))], [M.item((1,3*(k-1))), M.item((1,3*(k-1)+1)), M.item((1,3*(k-1)+2))], [M.item((2,3*(k-1))), M.item((2,3*(k-1)+1)), M.item((2,3*(k-1)+2))]])
        return uk
    else:
        return 'You entered a value of k out of range: k must be an integer between 1 and 3'

def comb(U1,U2,U3):
    """Combines three 3x3 matrices into a 3x9 matrix"""
    U=matrix([[U1[i,0], U1[i,1], U1[i,2], U2[i,0], U2[i,1], U2[i,2], U3[i,0], U3[i,1], U3[i,2]] for i in xrange(3)])
    return U

def numm(ss):
    """Transform triples of numbers into 3x9 matrices"""
    if isinstance(ss,matrix):
        return ss
    else:
        return comb(transfnm(ss[0]),transfnm(ss[1]),transfnm(ss[2]))

def matt(M):
    return [transfmn(extract(j,M%2)) for j in xrange(1,4)]

def p2(n):
    """Gives max{m:2^m<=n}"""
    j=0
    while 2**j<=n:
        j=j+1
    return j-1

def impp2(n): 
    """Gives [i_1,...,i_k]: n=2^(i_1)+...+2^(i_k)"""
    n1=[]    
    n1.extend([p2(n)])
    n2=n
    i=0
    while n2>2:        
        n2=n2-2**n1[i]
        if n2>0:
            n1.extend([p2(n2)])
        i=i+1        
    return n1

"""Default parameter of precision of inverse"""
precinv=1
#------------------------------------------------------------
#Define the class of infinite upper triangular matrices

class M:
    
    def __init__(self,k,*matrices):
        self.k = k
        matrices=[numm(i) % 2 for i in matrices]
        self.m = list(matrices)
        if not (not self.m):            
             while (self.m[0] == matrix ([[0for i in xrange(9)] for k in xrange(3)])). all():
                 """This guarantees that if the first matrices are zero matrices they are counted in the k and don't appear in the bracket"""             
                 self.m.pop(0) 
                 self.k = self.k +1
                 if not self.m:
                      self.k=0
                      break
        if not (not self.m):    
             while (self.m[len(self.m)-1] == matrix ([[0for i in xrange(9)] for k in xrange(3)])). all():
                 """This guarantees that if the last matrix is the zero matrix it does not appear in the argument"""
                 self.m.pop()
                 if not self.m:
                      self.k=0
                      break
        else:
            self.k=0
        
    def __str__(self):
        if not self.m:
            return 'Id'
        else:
            return 'M_%d (%s)' % (self.k,self.m)

    def __eq__(self,other):
        """A==B"""
        if len(self.m)==len(other.m)==0:
            return True
        else:
            for i in xrange (max(len(self.m),len(other.m))):
                if (len(self.m)!=len(other.m)) or (self.k!=other.k) or (self.m[i]!=other.m[i]).any():
                    return False
                    break
            return True

    def __ne__(self, other):
        """A!=B"""
        return not self.__eq__(other)

    def transfmn(self):
        """From matrix arguments to number arguments"""
        if not self.m:
            return 'Id'
        else:
            return 'M_%d (%s)' % (self.k, [[transfmn(extract(j,self.m[i])) for j in xrange(1,4)] for i in xrange(len(self.m)) ])

    def __repr__(self):
        """By just typing a name of an element g or an expression in the shell (without print) we get the equivalent of print g.transfmn() """
        return self.transfmn()
    
    def ext(self,n,m):
        """From an object in the class to the submatrix obtained by considering n rows and m columns"""
        SM=matrix ([[0 for i in xrange(m)] for k in xrange(n)],dtype=object)
        for i in xrange(n):
             for j in xrange(m):
                  if i==j:
                       SM[i,j]=matrix([[1,0,0],[0,1,0],[0,0,1]])                  
                  elif self.k<j-i<=len(self.m)+self.k: #Note: the nth upper diagonal is the locus of entries such that the the difference between the column and row is n
                       SM[i,j]=extract(i%3+1,self.m[j-i-self.k-1])
                  else:
                       SM[i,j]=matrix([[0,0,0],[0,0,0],[0,0,0]])                        
        return SM

    def __mul__(self,other):
        """A*B for A in M and B in M or MX"""
        if isinstance(other,M):
            """If B is in M"""
            res1=self.ext(3,self.k+len(self.m)+3)
            res2=other.ext(self.k+len(self.m)+3,other.k+len(other.m)+self.k+len(self.m)+4)
            res=(res1*res2)%2
            """Turn the result into the standard form M_d(....)"""
            res=[res.diagonal(i) for i in xrange(1,other.k+len(other.m)+self.k+len(self.m)+2)]
            res=[n.tolist() for n in res]
            A=[0 for i in xrange(len(res))]
            for i in xrange(len(res)):
                 A[i]=comb(*res[i][0])
            return M(0,*A)        
        else:
             """If B is in MX"""
             m=other.k+len(other.l)
             res1=self.ext(3,m+3)
             res2=other.ext(m+3,m+3)
             res=(res1*res2)%2
             """Turn the result into the standard form M_d(....)"""
             res=[res.diagonal(i) for i in xrange(1,m+1)]
             res=[n.tolist() for n in res]
             A=[0 for i in xrange(len(res))]
             for i in xrange(len(res)):
                A[i]=comb(*res[i][0])
             return MX(0,*A)

    def inv(self):
        """Matrix inverse: since it is not clear where to cat the matrix, we let parameter a range and break the loop when product gives Id"""
        for a in xrange(precinv*(self.k+len(self.m))+3, precinv*(self.k+len(self.m))+4):
            res1=self.ext(a,a)
            res=inversebigmatrix(res1,2)
            for n in xrange(a-1,2,-1):
                 res=np.delete(res,n,0)
            """Turn the result into the standard form M_d(....)"""
            res=[res.diagonal(i) for i in xrange(1,a-2)]
            res=[n.tolist() for n in res]
            A=[0 for i in xrange(len(res))]
            for i in xrange(len(res)):
                 A[i]=comb(*res[i][0])
            if (M(0,*A)*self==Id):# and (self*M(0,*A)==Id):
                 return M(0,*A)
                 break
        return MX(0,*A)
                
    def __pow__(self,other):
        """A^n where n is an integer-could be -ve"""
        if isinstance(self,M):
            if other>=0:
                K = M(0,)
                for i in range(other):
                   K = K * self
            else:
                K = M(0,)
                for i in range(-other):
                   K = K * (self.inv())
        return K

    def conj(self,other):
        """Given A,B returns A^(-1)*B*A"""
        if isinstance(other,M):
            return self**(-1)*other*self
        if isinstance(other,MX):          
            global precinv
            precinv0=precinv
            precinv=1+(other.k+len(other.l))/(self.k+len(self.m))
            Ainv=self**(-1)            
            precinv=precinv0
            return Ainv*other*self
    
    def comm0(self,other):
        """Commutators"""
        if isinstance(other,M):
            return (self.inv())*((other.inv())*(self*other))
        if isinstance(other,MX):
            global precinv
            precinv0=precinv
            precinv=1+(self.k+other.k+len(other.l)+1)/(self.k+len(self.m))
            Ainv=self**(-1)
            precinv=precinv0
            return Ainv*other.conj(self)

    def comm(self,other,*arg):
        if len(arg)==0:
            return self.comm0(other)
        else:
            t=list(arg)
            t.insert(0,self.comm0(other))
            while len(t)>=2:
                t[0]=(t[0]).comm0(t[1])
                t.pop(1)
            return t[0]

    def commr(self,other,*arg):
        if len(arg)==0:
            return self.comm0(other)
        else:
            t=list(arg)
            while len(t)>=2:
                t[len(t)-2]=(t[len(t)-2]).comm0(t[len(t)-1])
                t.pop()
            return self.comm0(other.comm0(*t))

    def __add__(self,other):
        """For A=M(k(A),a_1,...,a_m(A)), B=M(k(B),b_1,...,b_m(B)) returns A+B=M(min(k(A),k(B),....), similarly if B in MX"""
        if isinstance(other,M):
            a=max(self.k+len(self.m),other.k+len(other.m))
            res1=self.ext(3,a+3)
            res2=other.ext(3,a+3)
            res=res1+res2
            res=[res.diagonal(i) for i in xrange(1,a+1)]
            res=[n.tolist() for n in res]
            A=[0 for i in xrange(len(res))]
            for i in xrange(len(res)):
                A[i]=comb(*res[i][0])
            return M(0,*A)
        if isinstance(other,MX):
            a=other.k+len(other.l)
            res1=self.ext(3,a+3)
            res2=other.ext(3,a+3)
            res=res1+res2
            res=[res.diagonal(i) for i in xrange(1,a+1)]
            res=[n.tolist() for n in res]
            A=[0 for i in xrange(len(res))]
            for i in xrange(len(res)):
                A[i]=comb(*res[i][0])
            return MX(0,*A)

    def trunc(self,n):
        """Truncate an instance of M after n upperdiagonals-output is in MX"""
        if n<=self.k:
            return MX(n,)
        elif self.k<=n<=self.k+len(self.m):
            A=[self.m[i] for i in xrange(n-self.k)]
            return MX(self.k,*A)
        else:
            res=MX(self.k,*(self.m))
            for i in xrange (n-self.k-len(self.m)):
                (res.l).append(numm([0,0,0]))
            return res

        
    def truncf(self,funct,n,other,*arg):
        A=self.trunc(n)
        B=other.trunc(n)
        if not (not arg):
            arg1=[0 for i in xrange(len(arg))]
            for i in xrange(len(arg)):                
                arg1[i]=(arg[i]).trunc(n)
            return funct(A,B,*arg1)
        else:
            return funct(A,B)

def tr(f,n):
    def helper(*arg):
         arg1=[0 for i in xrange(len(arg))]
         for i in xrange(len(arg)):
             arg1[i]=(arg[i]).trunc(n)
         return f(*arg1)
    return helper


        
Id=M(0,)

#--------------------------------------------------
#Computations with matrices with unknown diagonals:

class MX:
    
    def __init__(self,k,*matrices):
        self.k = k
        matrices=[numm(i) % 2 for i in matrices]
        self.l = list(matrices)
        if not (not self.l):            
             while (self.l[0] == matrix ([[0for i in xrange(9)] for k in xrange(3)])). all():
                 """This guarantees that if the first matrices are zero matrices they are counted in the k and don't appear in the bracket"""
                 self.l.pop(0) 
                 self.k = self.k +1
                 if not self.l:                      
                      break
        
    def __str__(self):
        return 'M_%d (%s,?)' % (self.k,self.l)

    def __gt__(self,other):
        """A>B if A in MX, B in M or MX, B is strictly contained in A"""
        if isinstance(other,MX):
            if (not self.l):
                if self.k<other.k:
                    return True
                elif (self.k==other.k) and not(not other.l):
                    return True
                else:
                    return False
            else:
                if len(other.l)>len(self.l):
                    for i in xrange (len(self.l)):
                        if (self.k!=other.k) or (self.l[i]!=other.l[i]).any():
                            return False
                            break
                    return True
                else:
                    return False
        if isinstance(other,M):
            if (not self.l):
                if (self.k<=other.k) or (not other.m):
                    return True
                else:
                    return False
            else:
                if len(other.m)>=len(self.l):
                    for i in xrange (len(self.l)):
                        if (self.k!=other.k) or (self.l[i]!=other.m[i]).any():
                            return False
                            break
                    return True
                elif (len(other.m)<len(self.l)) and M(self.k,*(self.l))==other:
                    return True
                else:
                    return False

    def __ge__(self,other):
        """A>=B if A in MX, B in M or MX, B is contained or 'equal' to A"""
        if isinstance(other,MX):
            if self>other:
                return True
            elif (len(self.l)==len(other.l)==0) and (self.k==other.k):
                return True
            else:
                for i in xrange (max(len(self.l),len(other.l))):
                    if (len(self.l)!=len(other.l)) or (self.k!=other.k) or (self.l[i]!=other.l[i]).any():
                        return False
                        break               
                return True
        if isinstance(other,M):
            return self>other
        
    def __lt__(self,other):
        """A<B if A is in M or MX, B in MX, A strictly contained in B"""
        if isinstance(other,MX):
            if (not other.l):
                if other.k<self.k:
                    return True
                elif (self.k==other.k) and not(not self.l):
                    return True
                else:
                    return False
            else:
                if len(other.l)<len(self.l):
                    for i in xrange (len(other.l)):
                        if (self.k!=other.k) or (self.l[i]!=other.l[i]).any():
                            return False
                            break
                    return True
                else:
                    return False
        if isinstance(other,M):
            return False
        
    def __le__(self,other):
        """A<=B if A is in M or MX, B in MX, A is contained or 'equal' to B"""
        if isinstance (other,MX):
            if self<other:
                return  True
            elif (len(self.l)==len(other.l)==0) and (self.k==other.k):
                return True
            else:
                for i in xrange (max(len(self.l),len(other.l))):
                    if (len(self.l)!=len(other.l)) or (self.k!=other.k) or (self.l[i]!=other.l[i]).any():
                        return False
                        break               
                return True            
        if isinstance (other,M):
            return False
    
    def transfmn(self):
        """From matrix arguments to number arguments"""
        return 'M_%d (%s,?)' % (self.k, [[transfmn(extract(j,self.l[i])) for j in xrange(1,4)] for i in xrange(len(self.l)) ])

    def __repr__(self):
        """By just typing a name of an element g or an expression in the shell (without print) we get the equivalent of print g.transfmn() """
        return self.transfmn()

    def ext(self,n,m):
        """From an object in the class the submatrix obtained by considering n rows and m columns"""
        SM=matrix ([[0 for i in xrange(m)] for k in xrange(n)],dtype=object)
        for i in xrange(n):
             for j in xrange(m):
                  if i==j:
                       SM[i,j]=matrix([[1,0,0],[0,1,0],[0,0,1]])
                  elif self.k<j-i<=len(self.l)+self.k:
                       SM[i,j]=extract(i%3+1,self.l[j-i-self.k-1])
                  else:
                       SM[i,j]=matrix([[0,0,0],[0,0,0],[0,0,0]])        
        return SM

    def __mul__(self,other):
        """A*B for A in MX and B in M or MX"""
        if isinstance(other,M):
            """If B is is M"""
            m=self.k+len(self.l)
            res1=self.ext(3,m+3)
            res2=other.ext(m+3,m+3)
            res=(res1*res2)%2
            res=[res.diagonal(i) for i in xrange(1,m+1)]
        else:
            """If B is in MX"""
            if self==other:
                m1=self.k+len(self.l)
                m2=2*(self.k)+len(self.l)
                res1=self.ext(3,m1+3)
                res2=other.ext(m1+3,m2+4)
                res=(res1*res2)%2
                res=[res.diagonal(i) for i in xrange(1,m2+2)]
            else:            
                m=min(self.k+len(self.l),other.k+len(other.l))
                res1=self.ext(3,m+3)
                res2=other.ext(m+3,m+3)
                res=(res1*res2)%2
                res=[res.diagonal(i) for i in xrange(1,m+1)]
        res=[n.tolist() for n in res]
        A=[0 for i in xrange(len(res))]
        for i in xrange(len(res)):
             A[i]=comb(*res[i][0])
        return MX(0,*A)

    def inv(self):
        "Matrix inverse"""
        a=self.k+len(self.l)+3
        if a<=3:
            return MX(0,)
        else:
            res1=self.ext(a,a)
            res=inversebigmatrix(res1,2)            
            for n in xrange(a-1,2,-1):
                res=np.delete(res,n,0)
            res=[res.diagonal(i) for i in xrange(1,a-2)]
            res=[n.tolist() for n in res]
            A=[0 for i in xrange(len(res))]
            for i in xrange(len(res)):
                A[i]=comb(*res[i][0])            
            return MX(0,*A)           

    def __pow__(self,other):
        """A^n where n is an integer-could be -ve"""        
        if isinstance(self,MX):
            if other!=0:
                if other<0:
                    self=self.inv()
                    other=-other                    
                r=impp2(other)
                S=[]
                for i in xrange(len(r)):
                    K=self
                    for j in xrange(r[i]):
                        K=K*K
                    S.extend([K])
                K=S[0]
                for i in xrange(1,len(r)):
                    K=K*S[i]
            else:
                K=M(0,)
        return K
    
    def conj(self,other):
        """Given A,B returns A^(-1)*B*A"""
        if isinstance(other,MX):
            a=other.k+min(self.k+len(self.l)+1,len(other.l))
            A=M(self.k,*(self.l))
            B=M(other.k,*(other.l))
        if isinstance(other,M):
            a=other.k+self.k+len(self.l)+1
            A=M(self.k,*(self.l))
            B=other
        global precinv
        precinv0=precinv
        precinv=1+a/(self.k+len(self.l))
        Ainv=A**(-1)
        precinv=precinv0
        res=Ainv*(B*A)
        if isinstance(res,M):
            mat=[res.m[i] for i in xrange(min(a-res.k,len(res.m)))]
            if a-res.k>len(mat):
                for i in xrange(-(len(res.m)-a+res.k)):
                    mat.append(numm([0,0,0])) 
        else:
            mat=[res.l[i] for i in xrange(min(a-res.k,len(res.l)))]          
            if a-res.k>len(mat):
                for i in xrange(-(len(res.l)-a+res.k)):
                    mat.append(numm([0,0,0]))            
        return MX(res.k,*mat)


            
    def comm0(self,other):
        """Commutators"""
        if isinstance(other,MX):
            if self.k+len(self.l)<=other.k+len(other.l):
                A=MX(other.k,*(other.l))
                if len(other.l)<self.k+len(self.l)+1:
                    for i in xrange(min(self.k+1,self.k+len(self.l)+1-len(other.l))):
                        (A.l).append(numm([0,0,0]))
                return (self.conj(A**(-1)))*A
            else:
                A=MX(self.k,*(self.l))
                if len(self.l)<other.k+len(other.l)+1:
                    for i in xrange(min(other.k+1,other.k+len(other.l)+1-len(self.l))):
                        (A.l).append(numm([0,0,0]))
                return A**(-1)*(other.conj(A))
        if isinstance(other,M):
            global precinv
            precinv0=precinv
            precinv=1+(other.k+self.k+len(self.l)+1)/(other.k+len(other.m))
            Binv=other**(-1)
            precinv=precinv0
            return (self.conj(Binv))*other

    def comm(self,other,*arg):
        if len(arg)==0:
            return self.comm0(other)
        else:
            t=list(arg)
            t.insert(0,self.comm0(other))
            while len(t)>=2:
                t[0]=(t[0]).comm0(t[1])
                t.pop(1)
            return t[0]

    def commr(self,other,*arg):
        if len(arg)==0:
            return self.comm0(other)
        else:
            t=list(arg)
            while len(t)>=2:
                t[len(t)-2]=(t[len(t)-2]).comm0(t[len(t)-1])
                t.pop()
            return self.comm0(other.comm0(*t))

    def __add__(self,other):
        """For A=MX(k(A),a_1,...,a_l(A)), B=MX(k(B),b_1,...,b_l(B)) or B=M(k(B),b_1,...,b_m(B)) returns A+B=MX(min(k(A),k(B),....)"""
        if isinstance(other,MX):
            if self==other:
                return M(0,)
            else:
                a=min(self.k+len(self.l),other.k+len(other.l)) #Note that in the same fct in M there is max instead of min
                res1=self.ext(3,a+3)
                res2=other.ext(3,a+3)
                res=res1+res2
                res=[res.diagonal(i) for i in xrange(1,a+1)]
                res=[n.tolist() for n in res]
                A=[0 for i in xrange(len(res))]
                for i in xrange(len(res)):
                    A[i]=comb(*res[i][0])
                return MX(0,*A)
        if isinstance(other,M):
            return other+self

    def trunc(self,n):
        """Truncate an instance of MX after n upperdiagonals-output is in MX"""
        if n<=self.k:
            return MX(n,)
        elif self.k<=n<=self.k+len(self.l):
            A=[self.l[i] for i in xrange(n-self.k)]
            return MX(self.k,*A)
        else:
            return self        
        
        
"""Alpha_i, Beta_i, Gamma_i"""        
alpha1 = matrix( [[0,0,0,0,1,1,0,1,0],[0,1,0,1,0,0,0,0,1],[1,1,1,0,0,0,0,1,0]])
beta1 = matrix( [[0,0,0,0,0,1,0,1,1],[1,0,1,0,0,0,0,1,1],[1,1,0,1,0,0,0,0,1]] )
gamma1 = matrix ([[0,0,0,0,1,1,0,1,0],[0,1,1,1,0,1,0,0,0],[1,0,0,0,1,1,0,0,1]])
alpha2 = matrix( [[0,0,0,0,1,0,0,0,1],[0,0,1,0,1,0,1,0,0],[1,1,0,0,1,1,1,0,0]])
beta2 = matrix ( [[0,0,0,0,0,0,0,0,0],[0,1,1,0,1,1,0,1,1],[0,1,0,0,1,0,0,1,0]])
gamma2 = matrix ( [[0,0,0,0,1,1,0,1,0],[1,0,0,0,1,0,1,1,1],[1,1,1,0,0,0,0,1,0]])
alpha3 = matrix ( [[0,0,0,0,0,1,0,1,1],[0,0,0,1,0,1,1,1,0],[0,0,0,0,1,0,1,1,1]])
beta3 = matrix ( [[0,0,0,0,1,0,0,0,1],[0,0,0,0,1,1,1,0,1],[0,0,0,1,0,1,0,1,0]])

"""Generators"""
x0=M(0,[11,11,11],[17,17,17],[26,26,26], [11,11,0],[17,0,0])
x1=M(0,[23,224,138],[59,136,495],[26,488,227],[23,224,0],[59,0,0])
x2=M(0,[46,68,217],[12,194,363],[26,326,77],[46,68,0],[12,0,0])
x3=(x0*x1)**(-1)
x4=(x1*x2)**(-1)
x5=x0*x1*(x2**(-1))
x6=(x0*x2)**(-1)
b3=M(0,numm([28, 235, 129]), numm([46, 138, 451]), numm([11, 479, 197]),numm([28, 235, 0]), numm([46, 0, 0]))
b4=M(0,numm([57, 164, 83]), numm([50, 210, 497]), numm([11, 118, 418]), numm([57, 164, 0]), numm([50, 0, 0]))
b5=M(0,numm([50, 175, 88]), numm([23, 83, 400]), numm([11, 418, 184]), numm([50, 175, 0]), numm([23, 0, 0]))
b6=M(0,numm([37, 79, 210]), numm([57, 217, 329]), numm([11, 364, 118]), numm([37, 79, 0]), numm([57, 0, 0]))
              
x=MX(0,numm([28,235,129]),numm([29,211,263]))
y=MX(0,numm([28,235,129]),numm([58,3,445]))

#----------------------------------------------------
#Linear combinations of alphai,betai,gammai
def factors_set():
    factors_set = ( (i,j,k) for i in [0,1] 
                          for j in [0,1]
                          for k in [0,1])
    for factor in factors_set:
        yield factor

def factors_set2():
    factors_set2 = ((i,j) for i in [0,1]
                        for j in [0,1])
    for factor in factors_set2:
        yield factor

def lincomb(M):
    weighs1 = (alpha1,beta1,gamma1)
    weighs2 = (alpha2,beta2,gamma2)
    weighs3 = (alpha3,beta3)
      
    for factors in factors_set():
        sum = matrix([[0,0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0,0]])
        for i in range(len(factors)):
            sum += (factors[i] * weighs1[i]) %2
            sum=sum%2
        if (sum == numm(M) %2).all():
            we1=('alpha1','beta1','gamma1')
            res=''
            for i in xrange(3):                
                if factors[i]!=0:
                    res+='+%s' %(we1[i])
            if len(res)==0:
                return '0'
            else:
                return res[1:]
    for factors in factors_set():
        sum = matrix([[0,0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0,0]])
        for i in range(len(factors)):
            sum += (factors[i] * weighs2[i]) %2
            sum=sum%2
        if (sum == numm(M) %2).all():
            we2=('alpha2','beta2','gamma2')
            res=''
            for i in xrange(3):                
                if factors[i]!=0:
                    res+='+%s' %(we2[i])
            return res[1:]
    for factors in factors_set2():
        sum = matrix([[0,0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0,0]])
        for i in range(len(factors)):
            sum += (factors[i] * weighs3[i]) %2
            sum=sum%2
        if (sum == numm(M) %2).all():
            we3=('alpha3','beta3')
            res=''
            for i in xrange(2):                
                if factors[i]!=0:
                    res+='+%s' %(we3[i])
            return res[1:]

#Precision
commm0=M.comm

def trall(n):
    if n==0:
        M.comm=commm0
    else:
        M.comm=tr(MX.comm,n)